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ABSTRACT 

The source population of gamma-ray bursts (GRBs) declines towards the present epoch being consistent with 
the measured decline of the star formation rate. We show this using the brightness distribution of 3255 long 
BATSE GRBs found in an off-line scan of the BATSE continuous L024 s count rate records. The significance of 
this conclusion is enhanced by the detection of three GRBs with known redshifts brighter than 10^^ erg s"' during 
the last two years. This is an argument in favor of the generally believed idea that GRBs are strongly correlated 
with the star production, at least on cosmological time scales, and favors the association of long GRBs with 
collapses of supermassive stars. However, we still cannot rule out neutron star mergers if the typical delay time 
for binary system evolution is relatively short. If we assume a steep decline of the GRB population at z > 1 .5, then 
their luminosity function can be clearly outlined. The luminosity function is close to a power law, dN/dL (x L~' 
for low luminosities over at least 1 .7 orders of magnitude. Then the luminosity function breaks to a steeper slope 
or to an exponential decline around 3 ■ 10^' erg s"' in the 50 - 300 keV range assuming isotropic emission. 

Subject headings: gamma-rays: bursts - methods: data analysis 



I. INTRODUCTION 

In spite of the remarkable progress that has proved the cos- 
mological origin of GRBs, there remains a number of extremely 
important issues which are still not resolved. Some of the main 
issues dealt with in this work are the following, (i) What is the 
the cosmological evolution of the source population of gamma- 
ray bursts (GRBs)? Does it evolve as the star formation rate 
(hereafter the SF rate) or does it has its own specific evolution? 

(ii) What is the luminosity distribution (or function) of GRBs? 

(iii) What is the total rate of GRBs in the universe? 

There are two competing approaches for such studies: (i) the 
"statisitical" one using large samples of poorly localized GRBs, 
and (ii) the "individual" one, using the small sample of optically 
followed-up GRBs, where we have additional information per 
event including the redshift and the intrinsic luminosity. 

The main data array for statistical studies was supplied by 
the Burst And Transient Source Experiment (BATSE) (Fish- 
man et al. 1989) onboard the Compton Gamma-Ray Obser- 
vatory iCGRO). The BATSE sample"* is a few times larger than 
the yield of all other experiments that have detected GRBs. It 
includes 2702 events in its final form. 

Fitting the BATSE data to various cosmological/evolutionary 
models has been the subject of many studies since the start of 
the BATSE operation in 199L For references to early works 
on fitting the BATSE brightness distribution of GRBs to cos- 
mological models, see Bulik (1999). For a number of results 
setting an upper limit to the width of the luminosity function 
of GRBs, see, e.g., Hakkila et al. (1996). Then Loredo & 
Wasserman (1998) demonstrated that the 3rd BATSE catalog 
does not constrain the luminosity function. Later work using 
the larger sample of the 4th BATSE catalog gave few con- 
straints. Krumholtz, Thorsett, & Harrison (1998) fitting the 
BATSE sample at peak fluxes P > 0.42 photons s"' cm"^ found 
that both a non-evolving GRB source population and SF evo- 



lutionary models fit the data even using standard candle GRBs. 
Wijers et al. (1998) demonstrated the same for the SF model. 
Totani (1999) showed that the SF model does not fit the data 
using the standard candle assumption. Panchenko (1999) using 
the SF model accounting for the time delay due to binary sys- 
tem evolution estimated the minimal width of the luminosity 
function as being two orders of magnitude in luminosity. Fi- 
nally, the recent work of Porciany & Madau (2001) deals with 
a larger sample including the BATSE non-triggered bursts of 
Kommers et al. (2000) extending the peak flux down to 0.18 
photons s"' cm"^. They found that the data fit cannot distin- 
guish between different variants of the GRB source evolution 
at large redshifts. However, they did not check how sensitive 
the fit is to the evolution at low redshifts. 

The general impression arising from fitting the BATSE data 
using cosmological models was that this approach had little fu- 
ture. Indeed, the "statistical" approach demonstrated an agree- 
ment between the data and a wide set of models. The very few 
constraints obtained were trivial. 

The main reason for the poor progress so far is the insuffi- 
cient depth of the BATSE sample, i.e., the too narrow bright- 
ness range. The brightest burst has a peak photon flux of 160 
photons s"' cm"^. Bursts useful for a usual least fit ^"e^ how- 
ever, at P < 30 photons s"' cm"^. The BATSE trigger threshold 
is at 0.2 photons s"' cm"^ but the bursts near the threshold are 
difficult to use because of poorly known threshold effects. All 
works cited above used the peak flux range above 0.4 or even 
above 1 photon s"' cm"^, which is then narrower than two or- 
ders of magnitude. Loredo and Wasserman (1998) showed that 
this is a too narrow range for obtaining constraints from the fits. 

The "individual" approach gave a wealth of important data. 
By itself this approach is, however, still unable to resolve 
the issues stated above. Besides having a too poor statistics 
(presently we have only 17 GRBs with known redshifts), the 
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approach is subjected to very strong selection biases. Never- 
theless, this small sample tells us that the luminosity function 
is at least 2.5 orders of magnitude wide and that it extends up 
to ~ 3 • 10^^ erg s"^ 

After most previous work in this area was completed, the 
following progress concerning the data accumulation has taken 
place: 

- BATSE obtained additional data until the deorbiting of 
CGRO in June 2000. 

- Searches for non-triggered bursts were performed by 
Schmidt (1999), Kommers et al. (2000), and Stern et al. (2000, 
2001). In the latter work, the statistics of useful GRB events 
was increased by a factor 1 .7 and the threshold effect was mea- 
sured. Thus, the useful fitting range was extended down to 0. 1 
photons s~' cm"^. 

- A sample of GRBs with known intrinsic luminosities 
(presently 17 events) appeared due to the optical afterglow ob- 
servations. 

In this work, we use these advances including the "individ- 
ual" GRB data. In addition, we include the brightest peak flux 
interval, which is statistically poor, into the maximum UkeU- 
hood fit and find that this peak flux interval is very informative. 
Thus we extend the fitting brightness range to 3 orders of mag- 
nitude. This allow us to obtain a number of conclusive results. 

In §§2-4, we describe the set of fitted data, the fitting models 
including the cosmology, the source evolution, and the lumi- 
nosity function of GRBs, and finally the the fitting procedure. 
In §5, we present results of the fits and show that the scenario 
of a non-evolving population of GRBs does not fit the data. In- 
stead, we demonstrate that the GRB-population should decUne 
approximately as fast as the star formation rate. We also de- 
termine the approximate shape of the luminosity function and 
give an estimate for the lower limit of the total rate of GRBs as 
being 3000 GRBs per year in the visible universe (that is up to 
a reasonably large redshift). 

2. THE DATA 

Probing various cosmological and evolutionary models we 
fitted the sample of 3255 BATSE GRBs longer than 1 s found 
by Stern et al. (2000, 2001) in the off-hne scan of the 
BATSE continuous daily records in 1.024 s time resolution. 
This sample, which is selected from the catalog of Stern & 
Tikhomirova*", is essentially uniform and has a corresponding 
efficiency matrix (measured by a test burst method), which is 
needed when fitting the weak end of the log N - log P distri- 
bution. (Hereafter the term log N - log P distribution means 
the differential distribution of GRBs vs. the logarithm of the 
peak photon flux). We excluded short bursts (consisting of one 
1 .024 s bin) from the analysis for two reasons: (i) short and 
long bursts could be separate phenomena, and (ii) the sample 
is incomplete regarding short bursts as they have a lower detec- 
tion efficiency and a wrong brightness estimate in 1.024 s time 
resolution. By excluding one-bin events, we make our sample 
more homogeneous. 

The brightness distribution of the GRBs in this sample was 
fitted using a hypothetical brightness distribution folded with 
the detection efficiency matrix described in Stern et al. (2001). 
This matrix was obtained using a sample of ~ 1 1,000 artificial 

* see littp://www.astro. su.se/groups/head/grb_arcliive.html 
' see, e.g., http://www.aip.de/ jcg/grb.html 

* see http://www.asdc.asi.it/bepposax/ 

' http://ssl.berkeley.edu/ipn3/index.html 



test bursts, which were superimposed on the BATSE continuous 
records and then passed through the same procedure of search 
and processing as real GRBs. The efficiency matrix is approxi- 
mately given by 
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where is the expected and c„ the measured count rate in units 
of counts s~^ cm"^; E{Ce) = 0.70(1 -exp[-(Ce/Ceo)^])'^ is the ef- 
ficiency function with fitted parameters Cgo = 0.097 counts s"' 
cm"^, v = 2.34; the log-normal factor describes the relative error 
of the measured count rate, a = 0.09(0. 08/ce)'/^, and the selec- 
tion bias is crudely expressed as c^o = Ce-i-0.05exp(-Ce/0.05). 

In order to constrain the intrinsic luminosity function (here- 
after the luminosity function or the LP), we used the sample 
of gamma-ray bursts with measured redshifts^. We cannot in- 
fer the LF from this sample as it is subjected to strong selec- 
tion biases. This is demonstrated in section 3.3. The redshift 
data, however, give us a useful piece of information, i.e., the 
existence of very intrinsically bright GRBs. Three of the in- 
trinsically brigfliest bursts are: GRB990123, GRB991216, and 
GRB000131 (named by dates), with redshifts 1.6 (Djorgovski 
et al. 1999), 1.02 (Vreeswijk et al. 1999), and 4.5 (Andersen et 
al. 2000), respectively, and with BATSE peak fluxes 16.4, 67.5, 
and 6.3 photons s"' cm"^ in the 50 - 300 keV range, respectively 
(estimated using the BATSE catalog). If they were emitted at 
z = 1, their peak fluxes taking into account the "K-correction" 
(i.e., the correction due to the spectral redshift effects on a fixed 
spectral band of the detector) would be 45, 69, and 84 photons 
s"' cm"^, respectively, assuming the cosmological parameters 
(I^Mj^^a) = (1,0). Hereafter, we use the photon peak flux at 
redshift z=l, /, as a measure of the intrinsic brightness. We, 
furthermore, use the intrinsic brightness interval of these three 
events, / > 40 photons s"' cm~^, to constrain the LF when fit- 
ting the BATSE log N -logP distribution. The choice of the 
three brightest events for this purpose is somewhat arbitrary. 
We cannot use a much wider brightness interval because of a 
brightness dependent selection bias. On the other hand, we can 
neglect such problems for the narrow brightness range of these 
three events. The data we fitted are presented in Figure 1. 

In order to impose a proper constraint on the LF we must 
estimate the sampling function for strong GRBs. With the sam- 
pling function we mean the probability that a burst will be de- 
tected, locaUzed, its afterglow observed and its redshift mea- 
sured. This function evolves with time. It was zero before 
1997. Then this function was limited by the field-of-view of 
the two Beppo-Sax Wide Field Cameras, ~ 0.08 of the sky**, as 
this was the main instrument supplying precise coordinates of 
GRBs during 1997 - 1998. 

In 1999 and 2000, many precise localizations were made 
by other systems, with most of them being made by the in- 
terplanetary network UlysseslKonuslNear (see the IPN home 
page'). This means that in this period the sampUng function 
became larger and for very bright events it could, in principle, 
approach unity as all instruments of the IPN had a A-k field-of- 
view. Actually it should be considerably less as in the same 
period, three very strong BATSE events (triggers 7301, 7491, 
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7595) were not localized. For one strong BATSE event (trigger 
7954, GRBOOOl 15) an X-ray transient was found, but no opti- 
cal transient. With this background, let us take a conservatively 
high estimate of the sampling function, 540 = 0.5 for GRBs with 
/ > 40 photons s"' cm"^, and the conservatively low estimate 
of the rate of these GRBs, A'4o = 3 yr"', in the visible universe. 

3. FITTING MODELS 

3.1. Cosmology 

We tried two sets of cosmological parameters, the flat matter- 
dominated universe, commonly used in most of previous works: 
(Om, J^a) = (1,0) and the vacuum-dominated cosmology, which 
is supported by recent data (ilM,f^A) = (0.3, 0.7) (see, e.g., 
Lukash, 2000). Hereafter, these two models are denoted as 
M-models and A-models, respectively. The distribution of 
GRBs over redshift for a non-evolving (NE) population for 
Om + Oa = 1 is 



dN 



oc 



dz' 



The photon number "luminosity" distance is defined by 

dz' 



lo ^nA+nud+z')^ 
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where Hq is the Hubble constant, which is assumed to be 75 km 
s~' Mpc~\ when estimating the luminosities of GRBs. For the 
standard luminosity distance, the factor \/l+z is replaced with 
l+z. 

3.2. Evolution of the Source Population 

Another important component of the model is the evolution 
of the population of GRB sources. We probed four cases: a 
non-evolving population and three evolution functions corre- 
lated with the history of star formation following Porciani & 
Madau (2001). The declineing phase of the star formation (SF) 
rate at z < 1 .5 is a relatively well measured function of z. Its his- 
tory at z > 2 is, however, controversial. This issue is discussed 
in Porciani & Madau (2001) giving the relevant references. Be- 
low we reproduce three versions of the SF evolution suggested 
in that work: 



0.3e3'*^ _, 
''^SFi(z) = M© yr 'Mpc 



(e3.8^+45) 
i.e., decreasing SF at z > 1.5, 
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^SF2(z) = . Mq yr-^Mpc^^ 



(e3.4z + 22) 
i.e., roughly constant SF at z > 2, 
0.i34e3.05z 



^SF3(Z) = 
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i.e., increasing SF at large z. 

Models will from now on be denoted as NE,M; SF1,A; and so 
on, where M and A denotes the two types of cosmologies. 

The next step is the generation of the standard candle log A' 
- log P distributions. At this step, we introduced an additional 
broadening of the observed brightness distribution due to the 



K-correction which depends on the type of GRB spectrum. For 
this purpose, we obtained the log - log P distributions us- 
ing Monte-Carlo simulations. To each simulated GRB we pre- 
scribed one of 54 spectra of bright BATSE bursts parametrized 
by the Band expression (Band et al. 1993). All these template 
spectra were assumed to be emitted at z= 1. Then we sampled 
the z of the burst and the corresponding K-correction for the 50 
- 300 keV band was applied to the apparent brightness of the 
simulated GRB. 

The resulting log - log P distributions are shown in Figure 
2. If we neglect the K-correction which depends on the GRB 
spectrum, then each distribution in Figure 2 is a direct reflec- 
tion of the corresponding redshift distributions of GRBs. The 
standard candle luminosity (before applying the K-correction) 
corresponds to a peak photon flux of 1 photon s~^ cm~^ in the 
50-300keV bandatz=l. 

If we vary this value when fitting the data, none of the 8 mod- 
els (two cosmological cases, four evolutionary cases) would 
still fit the observed log - log P distribution of 3255 long 
BATSE GRBs (shown by the crosses in Figure 2). The fact that 
the evolution of the SFl type with the standard candle lumi- 
nosity function cannot fit data was shown by Totani (1999) and 
Lloyd & Petrosian (1999). The NE model gives the smallest 
deviation in this case, but still the value of is unacceptable 
(61 for 27 degrees of freedom). However, as it will be shown 
below, the discrepancy in the case of the NE model can not be 
compensated by any hypothesis of the luminosity function. 

The slopes of the three log A' - log P distributions for an 
evolving SF rate are close to the Euchdean slope -3/2 starting 
from a peak flux P ~ 1 photon s~^ cm~^, which at this standard 
candle brightness corresponds to z = 1 . The agreement with the 
Euclidean slope is an accidental coincidence due to the super- 
position of cosmological and evolutionary effects. 

3.3. Parametrization of the Luminosity Function 

As was stated above, the luminosity distribution of events 
with known z cannot be used as a base for the LF model in 
the whole brightness range. This fact is clear from Figure 3 
where we present log A' - log P distributions that would give 
the sample of 17 GRBs with known absolute luminosities for 
the NE, SF2 and SF3 evolutionary models. All models give a 
striking disagreement with the data. The only way to reduce 
this disagreement is to assume an unreasonably sharp increase 
of GRBs at large redshifts which, in turn, will contradict the 
redshift data. Therefore the shape of a hypothetical LF remains 
arbitrary (except the constraints imposed on the brightest end 
of the LF). 

In order to get a handle on the LF of GRBs, we tried dif- 
ferent types of functions that describe common shapes of wide 
distributions in nature: the log-normal distribution (LGN), a 
truncated power law (TPL), a power law with an exponential 
cutoff (PLexp), and a broken power law (BPL). 

LGN: dN/dl = C - exp(-(ln^(///o)/252), with fliree free pa- 
rameters S, and C. 

TPL: dN/dl = C ■ /""' for /, < / < k and outside this in- 
terval. Free parameters are a, h,Io, and C. In some fits, h was 
fixed to -00 leaving 3 free parameters. 

PLexp: dN/dl = C ■ /""^ • exp(-///o), with three free pareme- 
ters a, Iq, and C. 

BPL: dN/dl = C ■ for h < I < k, dN/dl = Ci ■ I"^^'^ 
for /o < / < /i and dN/dl = outside the [h,l2] interval. Free 
parameters are a, (3, h, Iq, and C, while h was fixed to a value 
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above the maximum observed GRB brightness. 

We also considered a smoothed version of the broken power 
law: 

SBPL: dN/dl = CI"-^ /( 1 + (lo/iy ) 

For technical convenience, we measure the intrinsic bright- 
ness as peak count rate or peak photon flux, /, in the 50 - 300 
keV range produced by a GRB at a distance corresponding to 
z = 1. The absolute peak luminosity of the GRBs is related to 
/ as L = / • 3 • lO^'' erg s~^ assuming isotropic emission. Below 
we present the main results on the LF both in I and in absolute 
luminosity units. 



4. THE FITTING PROCEDURE 

We used the forward folding method when fitting the ob- 
served distribution of GRBs, i.e., the hypothetical brightness 
distribution was convolved with the efficiency matrix (1) and 
fitted to the observed distribution of GRBs over peak count rate 
(crosses in Figure 2). This distribution was represented by 29 
data points spaced by 0. 1 in log P in the interval 0.067 - 50 pho- 
tons s~^ cm~^ in the 50 - 300 keV range. During 9.1 years, there 
are 3 GRBs detected by BATSE that are brighter (the rightmost 
cross in Figure 2 and one GRB at 160 photons s"'cm"^, which 
is not shown). We treat this brightness range separately, es- 
timating the likehhood function of the fit for each peak flux 
interval. For the main interval, this is the standard proba- 
bility function. For the tail of the brightness distribution, the 
likelihood is the Poisson probability to sample not more then 3 
events brighter than 50 photons s ' cm"^ at the given number 
of such events for the full observation period, M50, predicted by 
the model. The final likelihood finction is the product of these 
two factors. 

Instead of the usual minimization procedure, we explore 
the parameter space sampling ^10^ random points. This is suf- 
ficient to find the minimum of with a good accuracy while at 
the same time investigating the x^ "topography" of the param- 
eter space region. This method does not work if the "valley" of 
the minimum has a very small volume in parameter space. For 
our cases, the minima are wide and smooth enough. 

The maximum likelihood point for the whole sample of 
points in parameter space represents the unconstrained fit, 
where the requirement of the redshift data, A'40 > 3 yr"\ was ig- 
nored. The subsample selected including this requirement rep- 
resents the constrained fit. The effective number of degrees of 
freedom (DOF) in the second case is smaller by unity as com- 
pared to the unconstrained fit. We present results of both the 
constrained and the unconstrained fits in order to demonstrate 
the role of the redshift data and the possible effects of the uncer- 
tainty in the estimate of the rate of intrinsically strong GRBs. 

The best fit parameters are presented in Tables 1 and 2. Table 
1 includes fits using various LF models and selected cosmolog- 
ical models. Table 2 summarizes fits using a broken power law 
LF for all cosmological models and provides data to compare 
their relevance using a Bayesian approach. 

5. RESULTS 

5.1. Rejection of models without source evolution 

The best unconstrained fit using the NE models has x^ = 
36.8 for 25 degrees of freedom (Table 1, truncated power law, 
TPL), which is marginally acceptable if we ignore the tail of 
the brightness distribution. For the tail, the model predicts M50 



= 11.5, while the real number is 3 yr"^ and the correspond- 
ing maximum likelihood drops down to 1.3-10"^ (see Table 1). 
When we impose the redshift data constraint, A'40 > 3 yr"', the 
maximum likelihood factor that we can obtain using the NE,A 
model is 3.2- 10"*, which is for the case of a broken power law 
LF (see Table 2). For the NE,M model, the results are even 
worse. 

The Ukehhood factor is not yet a rejection factor for the NE 
model, because we cannot exclude some bias or contamination 
which would increase the x^. To reject the NE hypothesis, we 
should demonstrate a good fit for other equally simple and rea- 
sonable models. Indeed, data fits with the SF models are much 
better (see Tables 1 and 2). Their maximum likelihood factor 
is about 0.02. This is the case when we can apply the Bayesian 
approach. The estimate of the rejection level for models with 
nonevolving GRB source population is the ratio of the maxi- 
mum NE likehhood factor (NE,A model in Table 2) to that for 
SF models (e.g, SF1,A in Table 2), which is 1.4 • lO""*. 

Figures 4, 5, and 6 demonstrate the differences between fits 
using NE and SF models. As seen in Figure 5, the high bright- 
ness slope of the log A' - log P distribution for the broken power 
law NE model is too flat and it cannot be made considerably 
steeper by modifying the LF. Note that the NE standard candle 
log A' - log P distribution with 7=1 photon s"^ cm"^ is already 
flatter than the observed log A'^ - log P distribution. Then if the 
LF is extended to ~ 100 photons s~' cm"^, the model tail will 
be considerably flatter than the observed tail independently of 
how the LF is extended. 

Looking at the integral distribution of real GRBs in Figure 
6, one notices that it dechnes faster than an Euchdean distribu- 
tion. However, this is still not a statistically significant fact. The 
probability of such a deviation from the -3/2 slope by chance is 
0.1, which is relatively large (integral distributions are known 
to produce an iUusion of statistically significant features from 
fluctuations). Most probably we are just dealing with a moder- 
ate fluctuation. Nevertheless, the observed slope could really be 
steeper than the Euclidean one if the decline of GRBs is steep 
enough. More data are required to clarify this issue. 

The rejection of the NE models is significant even without 
the redshift data. The unconstrained rejection factor, taken 
as the ratio of the unconstrained likelihoods of the NE,A and 
SF1,A models is ~ 2 • 10"^ (see Table 2). This means that the 
result is not very sensitive to the value of the constraining A'40. 
If we, e.g., overestimate the rate of intrinsically strong GRBs 
by a factor 3 (suppose that it is a fluctuation) and that A'40 = 1, 
then the rejection factor is 4-10"'*. 

While the NE models are rejected using the low redshift be- 
havior, one does not obtain any preference for a certain kind of 
SF evolution at large redshifts. Different SF models give similar 
likelihood results (see the Tables) except for the SF3 scenario, 
which gives a slightly worse fit. Furthermore, the data do not 
allow us to distinguish between matter-dominated and vacuum- 
dominated cosmologies. 

5.2. The Shape of the Luminosity Function 

Once the NE models have been rejected at a significant level, 
we now concentrate on the SF models, i.e., models with evolu- 
tion of the GRB source population. There are two clear features 
of the LF that we see for all SF models: a near power law in- 
terval at the lower brightness range and a break or turnover to- 
wards the bright end of the distribution. Attempting to replace 
this construction of the LF by a log-normal LF gives a decrease 
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of the maximum likelihood by 2 orders of magnitude (see, e.g., 
the two SFl models in Table 1). 

If we study the topography for the broken power law LF, 
we find a power law fragment at least 1.7 orders of magnitudes 
wide and there is no upper limit on its width. It can be ar- 
bitrarily extended to lower brightnesses (see Fig. 7). Indeed, 
one can see from Figure 7 that the distribution reaches its 
asympthotics: by extending the LF further to lower brightnesses 
does not affect the model in the range of the data points. 

On the other hand, the slope of this power law is surprisingly 
well constrained, especially for the SFl model (Fig. 7), and it 
is sUgtly sensitive to the choosen SF behavior at large redshifts. 
The shape of the x^ minimum changes, however. Larger SF 
rates at high redshifts allow flatter slopes a (see Fig. 8). The 
X^ - q; plots for a matter-dominated cosmology, which are not 
shown, are very close to those shown in Figure 8. They are just 
slightly narrower and more symmetric. 

A break or a turnover is necessary at a high significance level. 
Its removal increases x^ by Ax^ ~ 30 for SFl and Ax^ ~ 19 
- 25 for SF2 models (see Table 2). The properties of the break 
are, however, less certain than the parameters of the power law 
fragment. The x^ minimum is non-parabolic, asymmetric and 
relatively wide as seen in Figure 9. All we can say is that some 
turnover in the power law LF is required at an intrinsic bright- 
ness of about 10 photons s~' cm"^ at z = 1, or 3 • 10^' erg 
s~' for isotropic emission. The fitted position of the break is 
slightly sensitive to the cosmological model but the difference 
is within the statistical errors. We, however, can not distinguish 
between a power law break and an exponential cutoff of the LF. 
Both give a good x^ and the same maximum UkeUhood (com- 
pare Tables 1 and 2). Summarizing, we can claim that the be- 
havior of the LF below / ~ 3 photons s~'cm~^ or 10^ 'erg s"' (in 
the 50 - 300 keV range) is close to the power law dN /dl oc /~'-^ 
(or, more exactly, the power law index a — I can vary from 
~-1.35 to ~-1.5 depending on the model). Then it breaks 
to a steeper slope. A smooth break (i.e., the parametrization 
SBPL in section 3.3) also gives a good result. The likelihood is 
0.016 for the SF2,A model when the difference in the slopes is 
large enough, i.e., /3 < -3. 

Figure 10 shows a set of the best fit LFs for different models. 
The reason why the best fit LFs is more or less well defined is 
clear when comparing the LF curves with the BATSE log A' - 
log P distribution shown by the crosses. The latter should be 
a convolution of the former with the standard candle curves in 
Figure 2. For the SFl and SF2 models, the standard candle 
curves are "narrower" then the BATSE log - log P distribu- 
tion. Therefore the LF should roughly correspond to the main 
features of the observed log A'^ - log P distribution: a power law 
with a turnover. In the case of the SFl model with narrower 
redshift distribution, the LF is closer to the observed brightness 
distribution (see Fig. 10). It is natural that the required turnover 
of the LF is sharper than that of the log N - log P distribution 
with the sharp LF turnover being smoothed by the convolution. 

5.3. The lower limit on the total GRB rate 

The last column of Table 2 shows the lower limits on the total 
rate of GRBs in the visible universe. These limits correspond 
to Ax^ = 7., i.e., to Icr of a distribution with 25 DOE. Note 
that their values are obtained using an abrupt cutoff of the LF 
at the dim end. For the SF2,A model, the highest allowed cut 
off is at ~ 0.4 • 10^" erg s~' . A more reaUstic smoother cut off 
would give a higher lower limit. We suggest that the value 3000 



GRBs per year gives a realistic estimate of the minimum GRB 
rate. This estimate is for long GRBs only. The result depends 
on the SF model in a natural way predicting a larger result for a 
higher SF rate at large redshifts. 

For the SF2,A model of the GRB source evolution, the high- 
est minimum comoving GRB rate at large redshifts is ~ 3 yr"' 
Gpc"^. At the present epoch, it becomes 0.13 yr~^ Gpc~^. This 
lower limit coincides with the value claimed by Porciani & 
Madau (2001) as the estimate derived from a fit of the data of 
Kommers et al. (2000). It is, however, very difficult to compare 
results as the slopes and the ranges of the LF are different. In 
fact, the shape of the LF derived by Porciany and Madau has 
the same parametrization as one of our models, a power law 
with an exponential cut off, but, neverheless, contradict our re- 
sults predicting a different slope and a finite estimate for the low 
brightness cutoff of the LF. The latter could be a consequence 
of the different behavior of the log A'^ - log P distribution near 
the threshold in Kommers et al. (2000). 

Our estimate is close to the result of Schmidt (2000) which 
is ^ 0.2 yr"' Gpc~^ at z = derived assuming that the present 
GRB rate is 10 times less than at z ~ 1 .5 (in our case the decUne 
is a factor 23). On the basis of this estimate, we can discuss 
the possible association of GRB980425 and the supernova SN 
1998bw (Bloom et al., 1999, see also Lamb, 1999 for a discus- 
sion). Schmidt (2000) estimates the probability of events such 
as GRB980425 if the association is real as being as low as 10"^ 
yr"'. Indeed, flie distance to SN 1998bw is ~ 40 Mpc. The 
corresponding sampling volume is 2.7-10"'' Gps-'. Our estimate 
of the GRBs rate 0.13 yr"' Gpc'^ with the cutoff ~ 3 • erg 
s"' . The 50 - 300 keV peak luminosity of GRB980425 is 3-10''^ 
erg s~' (our rough estimate). If the luminosity function extends 
down to this luminosity with the determined slope a- 1 « -1 .4, 
then the rate of GRBs above this luminosity will be ~ 20 times 
larger, i.e., 2.6 yr"' Gpc"^ at z = and a corresponding rate of 
0.7-10"^ yr~' in the samphng volume. Furthermore, one should 
take into account the probability that such events will be local- 
ized with the accuracy of a few arcseconds which reduces the 
estimate by at least an extra order of magnitude. Thus the prob- 
abUity of occurance and good localization of an event such as 
GRB980425 within ~ 40 Mpc with the power law behavior of 
LF down to 10"*^ erg s"' is below 10"^ per year Such a small 
probability can hardly be compensated by a break in the LF be- 
low the observational cut off. The break must be so steep that 
it will affect the observed log A^ - log P distribution near the 
BATSE flireshold. This is a strong argument that GRB980425 
might represent a different phenomenon that does not overlap 
with classic GRBs in its luminosity function (but at the same 
time being indistinguishable from typical GRBs in its general 
appearence). A more probable possibility is that we are dealing 
with an accidental coincidence. 

6. CONCLUSIONS 

The source population of GRBs sharply decUnes fromz ^ 1 .5 
towards the present epoch. This fact is established at the con- 
fidence level 10""*. The measured decline of the star formation 
rate being used as the evolution hypothesis for the source pop- 
ulation of GRBs fits the data satisfactorily. 

This is what is expected according to prevailing view of 
GRBs as being the product of stellar evolution. Such a scenario 
was already successfully applied for the description of the ob- 
served log A' - log P evolution, e.g., by Wijers et al. (1998), 
Panchenko (1999), Kommers et al. (2000), Schmidt (2000), 
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and Porciani & Madau (2001). However, nobody has quantita- 
tively demonstrated that an evolution similar to the SF rate at 
some range of z is necessary to describe data. We can mention 
just the work of Schmidt (2000) where the SF rate hypothesis 
fits the bright tail of the log A' - log P distribution better (by a 
visual impression) than the NE hypothesis. However, no quan- 
titative comparison of different scenarios has been done. We 
have here demonstrated that the decline of the GRB population 
consistent with the SF decline is a necessary requirement. 

The main issue is whether the GRBs are associated with the 
collapse of massive stars (coUapsars) or the merging of neutron 
stars (mergers). On the small scale, these scenarios differ re- 
garding the expected correlations of GRBs with star formation 
regions in galaxies: collapsar GRBs should be well correlated, 
merger GRBs should show no correlation with the star forma- 
tion. This fact stimulated searches for such correlations using 
optical GRB afterglows (see, e.g., the review of van Paradijs, 
Koveliotou, & Wijers, 2000). On the cosmological scale, such 
correlations should exist in both cases. However, in the merger 
scenario the occurance of GRBs will be described by the SF 
rate convolved with a delay function. If the latter extends to 
a few billions years, then the decline of the GRB population 
will be considerably flatter. For estimates of the delay func- 
tion for mergers, see Portegies Zwart & Yungelson (1998) and 
Panchenko et al. (1999). 

Qualitatively, our results favor the collapsar scenario. A very 
interesting situation occurs if the very steep, steeper than -3/2, 
slope of the log - log P tail as shown in Figure 5 will per- 
sist to larger P {Ulysses data reduced by Atteia, Boer & Hurley, 
1999 indicate that this can be the case). Then one will have 
to accept that the GRB progenitors are very massive collapsars 
whose population can decline faster than the general SF rate. 



But, at present, we do not have sufficient statistical argu- 
ments to rule out mergers. We believe that the data can give 
tight constraints on the delay function for the merger scenario 
and the latter can be challenged by such constraints. However, 
we leave such estimates for future studies for the reason that in 
order to obtain solid conclusions, it is worth to incorporate the 
Ulysses data, which would provide at least a doubling of the 
statistics of the brightest GRBs. 

Our results concerning the luminosity function of GRBs con- 
firm the conclusion of Loredo & Wasserman (1998) that the 
width of the luminosity function is not constrained by the 
BATSE data being wider than two orders of magnitude. The 
shape of the LF tried in the most of previous works was a trun- 
cated power law or a lognormal distribution which satisfied ear- 
lier data in a narrower brightness range. Schmidt (2000) used a 
broken power law hypothesis similar to one of our models. The 
position of the break is consistent with our results. 

The interpretation of the shape of the luminosity function is 
beyond the scope of this work. In principle, such type of distri- 
butions - broken or exponentially cut power laws are common 
in nature. Then the break implies some physical limit such as 
a finite energy source. We believe that the outlined shape of 
the LF might be a useful clue in the development of physical 
models for the GRB emission. 

This research made use of data obtained through the 
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work was supported by the Swedish Natural Science Research 
Council, the Royal Swedish Academy of Science, the Wenner- 
Gren Foundation for Scientific Research, and the Russian Foun- 
dation for Basic Research grant 00-02-16135. We thank the 
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Table 1 

Best hts of data to some cosmological models with different forms of the LF. 



Model 


LF 




M50 


Likelihood 


a 


S 


AIo 



NE,A 


TPL 


36.8,39.8 


11.5 


1.3-10^ 






2.8, 10.7 


ne,a 


LGN 


44.3,42.5 


16.6 


1.2-10"^ 




1.9 


0.032, 0.36 


ne,a 


PLexp 


54.0, 49.7 


18.7 


2-10"* 


-U.Zl_o.38 




15,21 


sfi,a 


LGN 


44.5,39.8 


6.0 


3- 10"^ 




2.0 


0.032, 0.12 


SFl.A 


PLexp 


31.2,30.7 


7.0 


0.022 


-0.43±0.09 




14,28 


SF2,A 


PLexp 


342,31.7 


6.7 


0.016 


-0.36_QiQ 




11,33 


SF2,M 


PLexp 


32.7,31.4 


6.9 


0.018 


-0 36**'^ 




12,26 



Note. — For models and LF parametrizations, see §§3.1-3.3. The fit with the LF unconstrained by 
the redshift data is given only for TPL case. For all other given fits, the LF is constrained by A'40 > 3 
yr"' . The values are given at the maximum likelihood and the minimum points, respectively. 
The "Likelihood" is the product of the x^ probability and the Poisson probability of sampling not 

more than 3 events at a given expected number M50. M50 is given at maximum Hkehhood point, a is 
the slope of the power law part of the LF. S is the width of the log-normal distribution (see §3.3). /o is 
the upper cut off brightness of the truncated power law LF, the center of the log-normal distribution, 
or the exponential cutoff energy for the PLexp LF (see §3.3). The errors in a and the confidence 
interval, A/q, for Iq correspond to Ax^ = 4 (corresponding to 217 in the conventional interpretation. 
However, one should be careful with such an interpretation when the results depends on missing data 
points below the threshold). 



Table 2 

Best hts of data to all evolution models with a broken power law LF 



Model x^ ^50 Likelihood Uncon. Ikh. a A/q Ax^ A'tot yr ^ 



NE,M 


39.1 


16.8 


6- 10'' 


4 • 10'' 


-0.32±0.50 


0.23, 3.3 


7. 


2800 


NE.,A 


40.4 


14.5 


3.2- 10~* 


0.87-10"^ 


"■•'-0,4 


0.43, 5.0 


5. 


3000 


SF1,M 


29.0 


7.2 


0.015 


0.043 


-0.49±0.08 


2.8, 9.2 


33 


1800 


SF1,A 


29.4 


6.7 


0.023 


0.050 


-0.48±0.09 


2.9, 16. 


28 


2000 


SF2,M 


30.2 


6.8 


0.018 


0.031 


-0.44±0.11 


2.9, 15. 


25 


2700 


SF2,A 


30.3 


6.0 


0.021 


0.028 


-0.43±0.12 


3.4,21. 


19 


3000 


SF3,M 


30.9 


7.2 


0.014 


0.20 


-0.38±0.12 


2.9,11. 


21 


4400 


SF3,A 


30.7 


6.55 


0.014 


0.17 




4.2, 21. 


16 


4500 



Note. — The LF is parametrized as a truncated broken power law with 4 free parameters (the BPL model in 
§3.3). The number of degrees of freedom is 25. Many details are given in the Note of Table 1. a is the slope of 
the power law part of the LF below the break. The errors correspond to a Ax^ = 4 confidence interval. A/o is 
the confidence interval (Ax^ = 4) of the break in the LF, /o, in photons s~' cm~^ at z = 1, Introducing a break in a 
single power law gives the x^ reduction, Ax^, in the next column. A'tot is the lower limit of the total GRB rate per 
year. 
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Fig. 1 . — Fitted data. The histogram is the raw count rate distribution of the sample of 3255 long GRBs found in the continuous BATSE records. The solid curve 

is the efficiency function E(Ce) in equation (1) measured using the test burst method. The discrepancy at the threshold results from the brightness bias (see eq. 1) 
and is removed using the efficiency matrix. The thick horizontal bar on the right tail of the distribution shows the count rate contribution of the three intrinsically 
brightest events with measured redshifts if they were emitting at z=l. The number of these events is renormalized to the 9.1 years of the BATSE observations. Their 
original number is 3 over ~ 2 years. 




0.01 




0.001 0.01 0.1 
Peak photon flux, 50- 



1 10 100 
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1000 



Fig. 2. — Standard candle log N - log P distributions for different models. The standard candle brightness corresponds to a photon peak flux of 1 photon s cm 
at z= 1. From top to bottom: no evolution (NE), SF3 (eq. 6), SF2 (eq. 5), SFi (eq. 4). SoUd curves - (QmiS^a) = (1,0), dotted curves - (Qm,^a) = (03,0.7). The 
crosses represent the observed log N - log P distribution of 3255 long BATSE GRBs. 
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Fig. 3. — Test for the luminosity distribution inferred from the sample of 17 eve nts with known redshifts and intrinsic luminosities. Log N - log P distributions 
obtained with the luminosity distribution of this sample are compared to the BATSE data (crosses) for different evolutionary models. Solid curve: NE model; dotted 
curve: SF2 model; dashed curve: SF3 model. The striking disagreement which can not be compensated with a reasonable redshift distribution shows that one can 
not rely on redshift data in the whole luminosity range because of a strong luminosity dependent selection bias. The used cosmology was (Qm, ^a) = (0.3, 0.7). 
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Fig. 4. — Fits of the bright tail of the log N - log P distribution using a broken power law LF. The value of x" versus the predicted number of bright events M50 
at > 3 yr"' . Actual value of M50 is 3. Parameters Ii,Io,a, and P are random. The models are NE,A (upper panel) and SF1,A (lower panel). 
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Peak photon flux, 50-300 keV, photons s"' cm"^ 



Fig. 5. — Best fits of data to models with no (upper panel, NE) and with SFl (lower panel) GRB source evolution. Crosses are observed data points corrected using 
the efficiency matrix (1), soUd curves - the models, dotted curves - the model LF (a broken power law), dashed line - the Euclidean -3/2 slope. The LF corresponds 
to a GRB luminosity distance at z = 1 with an arbitrary normalization of the rate. The cosmological model has (fiMjf^A) = (0.3,0.7). For fitting parameters, see 
Table 2. 




Fig. 6. — The same best fit functions as in Figure 5 but in integral form. Histogram - the raw peak count rate distribution in integral form (N(> P) - the number of 
GRBs with peak flux larger than P), solid curve - the NE,A model, dotted curve - the SF1,A model, and dashed line - the Euchdean distribution. The model curves 
are convolved with the efficiency matrix to correspond to the raw data in this figure. 
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Fig. 7. — The confidence area for parameters of the power law fragment of the broken power law LF. The model is SFl.A. The confidence area corresponds to 
Ax^ < 6.17, which formally corresponds to a 2 tr confidence interval. However, see Notes of Table 1. 
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Fig. 8. — The profile of the minima for the power law slope a. Left panels: for a broken power law LF; right panels: for the PLexp LF. From bottom to top: 

the SFl.A, SF2,A, and SF3,A models. 
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Fig. 9. — The characteristics of the power law break in the broken power law LF for the SF1,A model. Upper panel: versus the position of the break, Ig, in the 
intrinsic peak brightness scale. Lower panel: versus the power law slope difference. Other parameters are random. 
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Fig. 10. — Best fit luminosity functions for different models. Dotted curves: a broken power law and the PLexp for the SF2,A model; dashed curve: the broken 
power law for the SF2,M model; dash-dotted curve: the smoothly broken power law (SBPL, see section 3.3) for the SF2,A model; solid curves: a broken power law 
and the PLexp for the SF1,A model. The normalisation of the intrinsic brightness scale corresponds to the peak photon Hux being emitted at z= 1. Crosses show 
real data points versus the apparent brightaess, where the distance for each event is unknown. Thus the ordinate for real data has a different meaning. 



